function [tctab, vtab] = cursor_fits(varargin)
%
% fit selected data in the window to the specified function
% This routine assumes that you have recently, and on the current data set,
% used the cursors on the main datac screen to set a window (only the X
% axes count) for the fits. 
% The calling argument should be one of:
% 'singleexp' - one exp fit to data
% 'doubleexp' - two exponential decays fit to the data
% 'singleepsc' - single exp decay with 1-exp rise to power (second arg)
% 'doubleexpc' - double exp decay with 1-exp rise to power (second arg)
% June, 2008. Paul Manis
% 
global ALLCH DFILE
global CFIT_Result

CFIT_Result = []; % clear out the result first.
[xa, ya]=getxycursor(); %#ok<NASGU>
if(nargin > 0)
    functionid = varargin{1}; % Power with 1 expl decay, 2 is double
    if(nargin > 1)
        power = varargin{2};
    else
        power = 2;
    end;
    if(nargin > 2)
        channel = varargin{3};
    else
        channel = 1; 
    end;
else
    functionid = 'singleexp'; % single exponential fit via chebyshev
    power = 1;
    channel = 1;
end;
if channel == 1
    ochannel = 2;
else
    ochannel = 1;
end;

tctab=[];
vtab=[];
nrec=size(ALLCH{channel});

lt=find_times(xa);
tpnts = lt(1):lt(2);
if(any(lt == Inf))
    return;
end;
time = make_time(DFILE);
tbase=time(1,1:(lt(2)-lt(1)+1));
ntpnts = length(tbase);

% skip factor (for testing...) 0 implies every trace is done, 1 is every
% other, etc.
skipfac = 0; 
reclist = 1:(skipfac+1):nrec(1);
nr = length(reclist);
fiterr=zeros(nr,1);
cur = zeros(nr, ntpnts);
fitcur = zeros(nr, ntpnts);
goodfit  = ones(nr, 1);
vtab=zeros(nr,1);
t1=zeros(nr,1);
t2=zeros(nr,1);
t3=zeros(nr,1);
am=zeros(nr, 1);
% - note that below we seed successive fits with the results from the previous fit..
k=1;


% start by fitting an average of the traces
avg = mean(ALLCH{channel}(:,tpnts), 1);
fprintf(1, 'File: %s Records: [%d-%d]\n', DFILE.filename, DFILE.record(1), DFILE.record(end));
fprintf(1, 'Average:\n');
switch functionid
    case 'singleexp'
        laminit = [0, -1, 5]; % very initial guess
        laminit(1) = mean(mean(ALLCH{channel}(:,1:10)));
        [laminit, fiterr(k), fitcur(k,:), gf,  funcstr] = ...
            singleexp(tbase, avg, laminit, power);
    case 'doubleexp'
        laminit = [0, -1, 5, 20, 0]; % very initial guess
        laminit(1) = mean(mean(ALLCH{channel}(:,1:10)));
        [laminit, fiterr(k), fitcur(k,:), gf, funcstr] = ...
            doubleexp(tbase, avg, laminit, power);
    case 'singleepsc'
        laminit = [-1,    0.2,    3,    0]; % very initial guess
        laminit(4) = mean(mean(ALLCH{channel}(:,1:10)));
        [laminit, fiterr(k), fitcur(k,:),  gf, funcstr] = ...
            EPSC_fit1(tbase, avg, laminit, power);
    case 'doubleepsc'
        laminit = [-2,    0.2,    2,    25,   0,    -0.1]; % very initial guess
        laminit(6) = mean(mean(ALLCH{channel}(:,1:10)));
        [laminit, fiterr(k), fitcur(k,:),  gf, funcstr] = ...
            EPSC_fit2(tbase, avg, laminit, power);
    otherwise
        return;
end;
alam = zeros(length(reclist), length(laminit));

warning('off','MATLAB:divideByZero');

for j = reclist % note skip factor...
    i = nrec(1)-j+1;
    vstep = mean(ALLCH{ochannel}(i,tpnts(5:end-5))); %#ok<NASGU>
    QueMessage(sprintf('Fitting %d (%d of %d)', DFILE.record(i), i, nrec(1)));
    % reduce the number of points for a first pass fit
    if(k == 1)
        initpar = laminit; %#ok<NASGU>
    else % use previous fit for first guess.
        initpar = lam; %#ok<NASGU> % very initial guess
    end
    cur(k,:)=ALLCH{channel}(i,tpnts); % cut out only what we need
    vtab(k) = k;
    am(k) = max(cur(k,:));
    fprintf(1, 'Rec %d:  ', DFILE.record(i));

    switch functionid
        case 'singleexp'
 %           laminit(6) = mean(ALLCH{channel}(i, 1:10));
            initpar = laminit; % always seed from average
            [lam, fiterr(k), fitcur(k,:), goodfit(k)] = ...
                singleexp(tbase, cur(k,:), initpar, power);
            adc=lam(1); a0(k) = lam(2); t1(k)=lam(3); %#ok<AGROW,NASGU> % 
            %ts(k)=lam(6);
        case 'doubleexp'
 %           laminit(6) = mean(ALLCH{channel}(i, 1:10));
            initpar = laminit; % always seed from average
            [lam, fiterr(k), fitcur(k,:), goodfit(k)] = ...
                doubleexp(tbase, cur(k,:), initpar, power);
            adc=lam(1); a0(k) = lam(2); t1(k)=lam(3); t2(k)=lam(4); af1(k) = lam(5); %#ok<AGROW,NASGU>
            % ts(k)=lam(6);
        case 'singleepsc'
 %           laminit(6) = mean(ALLCH{channel}(i, 1:10));
            initpar = laminit; % always seed from average
            [lam, fiterr(k), fitcur(k,:), goodfit(k)] = ...
                EPSC_fit1(tbase, cur(k,:), initpar, power);
            adc=lam(4); a0(k) = lam(1); t1(k)=lam(2); t2(k)=lam(3); %#ok<AGROW,NASGU> 
            % ts(k)=lam(6);
        case 'doubleepsc'
            initpar = laminit;
            initpar(6) = mean(ALLCH{channel}(i, 1:10));

            [lam, fiterr(k), fitcur(k,:), goodfit(k)] = ...
                EPSC_fit2(tbase, cur(k,:), initpar, power);
            adc=lam(6); a0(k) = lam(1); t1(k)=lam(2); t2(k)=lam(3); 
            t3(k)=lam(4); af1(k)=lam(5); %#ok<AGROW,NASGU> % ts(k)=lam(6);
%            [f, fitcur(k,:)]=fit_func(lam, tbase, cur(k,:), 53, power, 0); % just to be sure at end

        otherwise
    end;
    alam(k,:) = lam;
    k = k + 1;
end
fprintf(1, 'Means: A0: %7.3f Tau(rise): %7.3f  Power: %d  Tau(decay): %7.3f\n',  ...
    mean(a0), mean(t1), power,  mean(t2))
CFIT_Result.functionid = functionid;
CFIT_Result.funcstr = funcstr;
CFIT_Result.adc = adc;
CFIT_Result.a0 = mean(a0);
CFIT_Result.taurise = mean(t1);
CFIT_Result.power = power;
CFIT_Result.tau1 = mean(t2);
if(strmatch(functionid, {'doubleexp', 'doubleepsc'}))
    fprintf(1, ' (double): Tau2(decay) %7.3f  Frac Tau2: %7.3f\n', mean(t3), mean(af1));
    CFIT_Result.tau2 = mean(t3);
    CFIT_Result.ftau2 = mean(af1);

end;

warning('on','MATLAB:divideByZero'); % restore warning status



%CONTROL(sf).FIT=FIT; % add fit to results.
%ctl = CONTROL;

h = findobj('Tag', 'Cursor_Fit'); % check for pre-existing window
if(isempty(h)) % if none, make one
    h = figure('Tag', 'Cursor_Fit', 'Name', 'Cursor Data Curve fits', ...
        'NumberTitle', 'off');
end
figure(h); % otherwise, select it
clf; % always clear the window...

orient landscape
% set(gcf, 'FontName', 'Arial');

fsize = 8;
msize = 3;
subplot('position', [0.15 0.32 0.30 0.20]);
plot(vtab, a0', '-ko', vtab, am, '-ks', 'MarkerSize', msize);
if(functionid == 2)
    hold on;
    plot(vtab, af1.*am', '-bs', vtab, (1-af1).*am', '-r^', ...
        'MarkerSize', msize);
end;
set(gca,'FontSize', fsize, 'FontName','Arial' );
%xlabel('V (mV)', 'Fontsize', 10);
ylabel('A (pA)', 'FontSize', fsize, 'FontName', 'Arial');

subplot('position', [0.15 0.075  0.30 0.20]);
plot(vtab, t1, '-ko', 'MarkerSize', msize);
set(gca, 'FontSize', fsize);
xlabel('V (mV)', 'FontSize', fsize);
ylabel('{\tau}{_r} (ms)', 'FontSize', fsize);

subplot('position', [0.55 0.075  0.30 0.20]);
plot(vtab, t2, '-gx', 'MarkerSize', msize);
set(gca, 'FontSize', fsize);
xlabel('V (mV)', 'FontSize', fsize);
ylabel('{\tau}{_f} (ms)', 'FontSize', fsize);

if(strcmp(functionid,{'doubleexp', 'doubleepsc'}))
    hold on;
    plot(vtab, t3, '-cx', 'MarkerSize', msize);
end;


subplot('position', [0.55 0.32  0.30 0.20]);
plot(vtab, fiterr, '-gv', 'MarkerSize', msize);
set(gca, 'FontSize', fsize);
ylabel('Error (pA)', 'FontSize', fsize);

subplot('position', [0.15, 0.57, 0.7, 0.35]);
set(gca, 'FontSize', fsize);

for i = 1: size(fitcur, 1)
    if(goodfit(i))
        fitcolor = 'green';
        trmark = '.';
        trcolor = 'black';
        trsize = 1;
    else
        fitcolor = 'red';
        trmark = 'o';
        trcolor = 'red';
        trsize = 5;
    end;
line(tbase, cur(i,:),...
    'color', trcolor,...
    'Marker', trmark,...
    'MarkerSize', trsize,...
    'LineStyle', 'none');

line(tbase, fitcur(i,:),...
    'color', fitcolor);
end;


title(sprintf('File: %s  Recs: [%d-%d] Window: %7.1f-%7.1f\n%s fits', ...
    DFILE.filename, DFILE.frec, DFILE.lrec,  xa(1), xa(2), functionid ),'FontSize',7);

return


function [lam, fiterr, fitcur, goodfit, funcstr] = EPSC_fit2(x_data, y_data, ...
    lam, alpha)
% Fit an EPSC to rise^power*(exp1 + exp2)

model=53;
funcstr = 'f.adc+f.a0.*((1-exp(-t/f.taurise)).^f.power) .* ((1-f.ftau2).*exp(-t/f.tau1)+f.tau2.*exp(-t/f.tau2))';
FitData(:,1) = x_data';
FitData(:,2) = y_data';

% pars:    A0      taurise   tau1  tau2 f_tau2    DC
pmask =   [ 1,         1,     1,        1,    1,    0];
lbound =  [-10000,      0.05,   0.5,        5,    0,   -2000];
ubound =  [ 10000,        10,     7,      200,    1,    2000];
order=length(lam)-1;
maxiter = 1500;
beta = 0;

watchon;
warning('off','MATLAB:divideByZero');
warning('off','MATLAB:nearlySingularMatrix');
[c,lam,goodfit]=curve_fitting(FitData(:,1), FitData(:,2), ...
    'levenberg','cubic', model, order, lam, pmask, lbound, ubound, alpha, beta, maxiter);
watchoff;
warning('on','MATLAB:divideByZero');
warning('on','MATLAB:nearlySingularMatrix');

[f, fitcur]=fit_func(lam, x_data, y_data, model, alpha, 0);
fiterr = norm(f);

result = ...
    sprintf('A=%8.3f Tr=%7.3f Tf=%7.3f Ts=%7.3f FracSlow=%7.3f  DC=%7.2f err=%4.3f',...
    lam(1), lam(2), lam(3), lam(4), lam(5), lam(6), fiterr);
fprintf(1, '%s\n', result);
return;


%------------------------- single exponential decay on epsc

function [lam, fiterr, fitcur, goodfit, funcstr] = EPSC_fit1(x_data, y_data, ...
    lam, alpha)
% Fit an EPSC to rise^power*(exp1 + exp2)

model=52;
funcstr = 'f.adc+f.a0.*((1-exp(-t/f.taurise) .* exp(-t/f.tau1))';

FitData(:,1) = x_data';
FitData(:,2) = y_data';

% pars:    A0      taurise       tau1   DC
pmask =   [ 1,    1,      1,     1];
lbound =  [-5000,  0.1,      1,    -2000];
ubound =  [ 5000,   10,     10,     2000];
order=length(lam)-1;
maxiter = 1500;
beta = 0;

watchon;

[c,lam, goodfit]=curve_fitting(FitData(:,1), FitData(:,2), ...
    'levenberg','cubic', model, order, lam, pmask, lbound, ubound, alpha, beta, maxiter);

watchoff;

[f, fitcur]=fit_func(lam, x_data, y_data, model, alpha, 0);
fiterr = norm(f);

result = sprintf('A=%8.3f Tr=%7.3f Tf=%7.3f  DC=%7.2f err=%4.3f',...
    lam(1), lam(2), lam(3), lam(4),  fiterr);
fprintf(1, '%s\n', result);
return;

%------------------------- single exponential decay

function [lam, fiterr, fitcur, goodfit funcstr] = singleexp(x_data, y_data, ...
    lam, alpha)
% Fit an single exponential decay to the data

model=50;
funcstr = 'f.adc+f.a0.*(exp(-t/f.tau1) ';

FitData(:,1) = x_data';
FitData(:,2) = y_data';

% pars:    DC     A0   Tau1
pmask =   [1,     1,      1  ];
lbound =  [-10000,  -20000     0.05];
ubound =  [10000,    20000,     200];
order=length(lam)-1;
maxiter = 1500;
beta = 0;

watchon;

[c,lam, goodfit]=curve_fitting(FitData(:,1), FitData(:,2), ...
    'levenberg','cubic', model, order, lam, pmask, lbound, ubound, alpha, beta, maxiter);

watchoff;

[f, fitcur]=fit_func(lam, x_data, y_data, model, alpha, 0);
fiterr = norm(f);

result = sprintf('DC=%8.3f A0=%7.3f  tau1=%7.2f err=%4.3f',...
    lam(1), lam(2), lam(3),  fiterr);
fprintf(1, '%s\n', result);
return;

%------------------------- double exponential decay

function [lam, fiterr, fitcur, goodfit, funcstr] = doubleexp(x_data, y_data, ...
    lam, alpha)
% Fit double exp decay to the data

model=51;


funcstr = 'f.adc+f.a0.*(1-f_tau2).*(exp(-t/f.tau1) .* f_tau2.*exp(-t/f.tau2))';

FitData(:,1) = x_data';
FitData(:,2) = y_data';

% pars:    DC A0      tau1    tau2   f_tau2 
pmask =   [1,     1,      1,      1,     1];
lbound =  [-10000,  -20000, 0.05,     5,     0];
ubound =  [10000,    20000,   20,     200,   1];
order=length(lam)-1;
maxiter = 1500;
beta = 0;

watchon;

[c,lam, goodfit]=curve_fitting(FitData(:,1), FitData(:,2), ...
    'levenberg','cubic', model, order, lam, pmask, lbound, ubound, ...
    alpha, beta, maxiter);

watchoff;

[f, fitcur]=fit_func(lam, x_data, y_data, model, alpha, 0);
fiterr = norm(f);

result = sprintf( ...
    'DC=%8.3f A0=%7.3f tau1=%7.3f  tau2 = %7.3f  fractau2=%7.3f err=%4.3f',...
    lam(1), lam(2), lam(3), lam(4), lam(5), fiterr);
fprintf(1, '%s\n', result);
return;


